clear all; close all; clc; 
% function compareStates
cucd=pwd; 

%% Paths to Load and Associated Names 
loadPath.A={'BMJ\DtrFG\FG4 oneOnly NoCorrel 89\End 07Q3 11 01 13'};
models.names.A='CperH'; 

loadPath.B={'BMJ\DtrFG\FG4 oneOnly NoCorrel AHE 89\End 07Q3 11 20 13'}; 
models.names.B='AHE'; 

loadPath.C={'BMJ\DtrFG\FG4 both NoCorrel\End 07Q4 12 03 13'}; 
models.names.C='Both Mode 1'; 

loadPath.D={'BMJ\DtrFG\FG4 both NoCorrel\End 07Q4 2ndMode 12 03 13'}; 
models.names.D='Both Mode 2'; 

%===========================================================================
%% What to Plot and Decompose 

%% 1. *stateDecomp* States to plot and decompose
compareStru.stateDecomp={'real rate fpfw','real rate gap','gdp gap'}; 
compareStru.stateDecompFig={'Real Rate FPFW','Real Rate Gap','GDP gap'}; 
% Number of variance decompositions 
compareStru.NVarDec=20; 

%% 2. *stateIRF* Impulse responses 
compareStru.stateIRF={'gdp lev','cons lev','invest lev','wages lev','hours','mgcost','inflation',...
    'ffr','real rate fpfw','real rate gap','gdp gap','wage gap','ffr lead 1','ffr lead 4'};
compareStru.stateIRFFig=compareStru.stateIRF; 
compareStru.NIRF=24; 

%% 3. *stateSpectrum* 
compareStru.stateSpectrum={'gdp lev','cons lev','invest lev','wages lev','hours','mgcost','inflation',...
    'ffr','real rate fpfw','real rate gap','gdp gap','wage gap','ffr lead 1','ffr lead 4'}; 
    %'Compensation per Hour','Average Weekly Earnings'};

%% 4. *stateMom* 
compareStru.stateMom={'gdp','cons','invest','hours','mgcost','inflation',...
    'ffr','real rate fpfw','real rate gap','gdp gap','wage gap','ffr lead 1','ffr lead 4'}; 
%'Compensation per Hour','Average Weekly Earnings'};
compareStru.stateMomFig=compareStru.stateMom; 
compareStru.NRep=1000; 
compareStru.NCorrel=4;
% if == nan will be set to [100 size(Y,1)]; 
compareStru.NVecSim=[100 80]; 

savePath=cr_dir( cucd,'BMJ','DtrFG',['Compare ',strdate]);     

%===========================================================================


%% Assign default values for the moments 
[~,compareStru]=ch_field(compareStru,'NSimVec',[100 150]);
[~,compareStru]=ch_field(compareStru,'NRep',100);
[~,compareStru]=ch_field(compareStru,'NCorrel',4);
[~,compareStru]=ch_field(compareStru,'NIRF',8);
[~,compareStru]=ch_field(compareStru,'NVarDec',0);
[~,compareStru]=ch_field(compareStru,'flagUnitIRF',0);

loadPath=orderfields(loadPath); 
compareStru.fields=fieldnames(loadPath); 

%% Dimensions, part I 
dim.models=length(compareStru.fields); 
dim.states=length(compareStru.stateDecomp); 
dim.IRFs=length(compareStru.stateIRF); 
dim.spectrum=length(compareStru.stateSpectrum);

compareStru.modelsCell=cell(dim.models,1);               
%==========================================================================
for ii=1:dim.models 
    
    %% Set *fldnm*, names and load paths 
    fldnm=char(compareStru.fields{ii});         
    
    disp(['Loading from :',char(loadPath.(fldnm))] ); 
    
    cd( fullfile(cucd,char(loadPath.(fldnm)) ) ); 
        
    disp(['Name: ',char(models.names.(fldnm))]); 
    compareStru.modelsCell(ii)={char(models.names.(fldnm))}; 
    
    
    %% Load model structures, assign output to models.<name>
    load modelStructures
    
    cd(cucd); 
    
    models.(fldnm).KFStru=KFStru; 
    models.(fldnm).names=names; 
    models.(fldnm).sampleStru=sampleStru; 
    models.(fldnm).report=report; 
    models.(fldnm).model=model; 
    
    if ii==1 
        %% Dimensions Part II 
        dim.y=size(model.ZZ,1);
        % Full state of transition equation
        dim.s=size(model.GG,1);
        % Innovations
        dim.eta=size(model.cholV,1);
        % Number of observations 
        dim.obs=length(sampleStru.vector);         
        % Assign sample vector 
        sample.vector=models.(fldnm).sampleStru.vector;
        % dimension for simulation  
        if isnan(compareStru.NVecSim)==true
            compareStru.NVecSim=[100 size(KFStru.smoothSt,1)];
            dispaj('Simulation size set to ',compareStru.NVecSim(1),' and ',compareStru.NVecSim(2));
        end
        
    else 
        if isequal( sample.vector , models.(fldnm).sampleStru.vector )==false 
            error('Samples do not match') 
        end 
    end 
    
    %% states.mat: Matrix with the states to compare 
    pos.extract=cellposition(compareStru.stateDecomp,models.(fldnm).names.states); 
    states.mat(:,:,ii)=models.(fldnm).KFStru.smoothSt(:,pos.extract);     
    
    clear model names report sampleStru KFStru; 
    
    compareStruTemp=compareStru; 
    
    %% reportPositions: assigns dimensions for what must be reported 
    compareStruTemp=reportPositions(compareStru,models.(fldnm).names,models.(fldnm).sampleStru);  
%      models.(fldnm).reportPos=reportPos; 
%      models.(fldnm).reportFlags=reportFlag; 
%     % clear reportPos reportFlag

     %      Not needed anymore, output included in reportPositions.m 
     %      momentStru=compareStru;
     %      momentStru=momentsAssignDefaults(momentStru,report,nameStru);

     %% Compute All moments 
     momentOut=momentsDraws(models.(fldnm).model.GG(:,:,end),...
         models.(fldnm).model.CC(:,end),...
         models.(fldnm).model.RR(:,:,end),....
         models.(fldnm).model.cholV(:,:,end),...
         models.(fldnm).model.ZZ(:,:,end),compareStruTemp);
     
     %% Table with Parameters and fit 
     
end 
cd(cucd); 

figStru=[];
figStru=plotAssignDefaults(figStru); 
figStru.legend=compareStru.modelsCell;
figStru.xaxis =sample.vector;
states.figHandle=zeros(dim.states,1); 
for ii=1:dim.states; 
    states.figHandle(ii)=figure; 
    tempMat=squeeze( states.mat(:,ii,:) ); 
    plot3samp(figStru.xaxis,tempMat,figStru); 
    title(char( compareStru.stateDecomp{ii} ) ); 
end 
    
cd(savePath); 
for jj=1:dim.states 
print(states.figHandle(jj),'-dpsc','statesFig','-append');
end 
cd(cucd); 
disp('Done'); 
    
    
          


    


